#-----------------------------------------------------------------------------#
#### FUNCTIONS ####
#-----------------------------------------------------------------------------#
gc()
library("lfe")
library(texreg)
library(R2HTML)
library(data.table)
library(readstata13)
library(Hmisc)
library(haven)
library(plyr,include.only = "rbind.fill")
memory.limit(64000)

simplelag <- function (x,by=NULL,mylag=1,outside=NA) {
  myend <- length(x)- mylag
  if(!is.null(by))  {
    lby <- c(replicate(mylag,""),as.character(by[1:myend]))
    y0 <- c(replicate(mylag,outside),x[1:myend])
    y <- ifelse(as.character(by)==lby,y0,outside)
  }
  else {
    y <- c(replicate(mylag,outside),x[1:myend])
  }
}



forward <- function (x,by=NULL,myforward=1,outside=NA) {
  mystart <- myforward+1
  if(!is.null(by))  {
    fby <- c(as.character(by[mystart:length(x)],replicate(myforward,"")))
    y0 <- c(x[mystart:length(x)],replicate(myforward,outside))
    y <- ifelse(as.character(by)==fby,y0,outside)
  }
  else {
    y <- c(x[mystart:length(x)],replicate(myforward,outside))
  }
}


retain <- function(x,event,outside=NA)
{
  indices <- c(1,which(event==TRUE),length(x)+1)
  values <- c(outside,x[event %in% TRUE])
  y <- rep(values,diff(indices))
}

obtain <- function(x,event,outside=NA)
{
  indices <- c(0,which(event==TRUE),length(x))
  values <- c(x[event %in% TRUE],outside)
  y <- rep(values,diff(indices))
}



wtd.summary <- function (x,w=NA) {
  name <- deparse(substitute(x))
  nb_obs <- length(x)
  sum_wgt <- sum(w[is.na(x) ==F])
  nb_missing <- sum(is.na(x))
  sum_wgt_missing <- sum(w[is.na(x) ==T])
  wtd_mean <- wtd.mean(x,w=w,na.rm=T)
  wtd_sd <- wtd.var(x,w=w,na.rm=T)**0.5
  min <- min(x,na.rm=T)
  p1 <- wtd.quantile(x,weights=w,probs=0.01)
  p10 <- wtd.quantile(x,weights=w,probs=0.10)
  q1 <- wtd.quantile(x,weights=w,probs=0.25)
  q2 <- wtd.quantile(x,weights=w,probs=0.50)
  q3 <- wtd.quantile(x,weights=w,probs=0.75)
  p90 <- wtd.quantile(x,weights=w,probs=0.90)
  p99 <- wtd.quantile(x,weights=w,probs=0.99)
  max <- max(x,na.rm=T)
  structure(data.frame(name,nb_obs,sum_wgt,nb_missing,sum_wgt_missing,wtd_mean,wtd_sd,min,p1,p10,q1,q2,q3,p90,p99,max))
}


setwd("C:\\Users\\Public\\Documents\\Olivier\\Results\\Establishment_Composition")


#-----------------------------------------------------------------------------#
#### DATA CREATION ####
#-----------------------------------------------------------------------------#

# # CONSTITUTION PROGRAM
# # this is a database with only exits and layoffs
# ll <- read_sas("../Layoffs/layoffs.sas7bdat",col_select = c("siret","year","pds","nblayoffs","nbexits","effref"))
# ll <-ll[ as.numeric(ll$year)>2001,]
# 
# # # etab is a database with all establishment surveyed in MMO/DDMO
# # etab <- data.frame(read_sas("../Layoffs/etab_mmo.sas7bdat",col_select=c("siret","pds","year","effref")))
# # etab <- etab [as.numeric(etab$year)>2001,]
# # saveRDS(etab,"../Layoffs/etab2002.RDS")
# etab <- readRDS("../Layoffs/etab2002.RDS")
# 
# ll <- merge(ll,etab,by=c("siret","year"),all=T)
# 
# ll$pds <- ifelse(is.na(ll$pds.x + ll$pds.y)==F, (ll$pds.x + ll$pds.y)/2,
#                  ifelse(is.na(ll$pds.x)==F, ll$pds.x,
#                         ifelse(is.na(ll$pds.y)==F, ll$pds.y,NA)))
# 
# ll$effref <- ifelse(is.na(ll$effref.x + ll$effref.y)==F, (ll$effref.x + ll$effref.y)/2,
#                  ifelse(is.na(ll$effref.x)==F, ll$effref.x,
#                         ifelse(is.na(ll$effref.y)==F, ll$effref.y,NA)))
# 
# ll$nbexits[is.na(ll$nbexits)] <- 0
# ll$nblayoffs[is.na(ll$nblayoffs)] <- 0
# str(ll)
# 
# ll <- ll[,-which(names(ll) %in% c("pds.x","pds.y","effref.x","effref.y"))]
# str(ll)

# View(ll[is.na(ll$effref),])
# sum(ll$nblayoffs)/sum(ll$nbexits)
# saveRDS(ll,"../Layoffs/layoffs2002.RDS")
ll <- readRDS("../Layoffs/layoffs2002.RDS")
# rm(ll)



ee<-readRDS("est.rds")
ee <- ee[ee$year>2001 & ee$year<2018,]
         
         
ll$r_lyof <- ifelse(ll$effref>0,ll$nblayoffs/ll$effref,NA)
# ll$r_otlyof <- ifelse(ll$effref>0,ll$nbothlayoffs/ll$effref,NA)
ll$r_exit <- ifelse(ll$effref>0,ll$nbexits/ll$effref,NA)
ll$year <- as.numeric(ll$year)

str(ee)
gg <- merge(ee,ll[,c("siret","year","nblayoffs","nbexits","effref","pds")],
            by.x=c("est","year"),by.y=c("siret","year"))

str(ll)
str(gg)

 
# gg$nblayoffs[gg$year %in% 2002] <- 0
# gg$nbexits[gg$year %in% 2002] <- 0
# 
# gg$r_lyof[gg$year %in% 2002] <- 0
# gg$r_exit[gg$year %in% 2002] <-0
 
table(gg$year)
str(gg)

rm(ff,ff2,ff3,e02,l02)
rm(aa,dd,ee,ll,nn,etab)

gc()

gg$f9010 <- gg$f9910+gg$f9099


gg$f9010xf9010 <- ifelse(gg$f9010>0,(gg$f9010-1)/(gg$count-1),NA)
gg$f9910xf9910 <- ifelse(gg$f9910>0,(gg$f9910-1)/(gg$count-1),NA)

gg_b <- gg[,c("est","year","count")]
gg_b$year <- gg_b$year + 1
colnames(gg_b)<- c("est","year","lcount")

gg <- merge(gg,gg_b,by=c("est","year"),all.x=T)
rm(gg_b)
gg <- gg[order(gg$est,gg$year),]

gg$lndnbwkrs_neg <- (log(gg$count)-log(gg$lcount))*((log(gg$count)-log(gg$lcount))<=0)
gg$lndnbwkrs_neg <- ifelse(is.na(gg$lndnbwkrs_neg),log(gg$count),gg$lndnbwkrs_neg)
gg$lnnbwkrs_cumneg <- ave(gg$lndnbwkrs_neg,gg$est,FUN=function(x) cumsum(x)) 

gg$sf9010 <- ave(gg$f9010,gg$year,FUN=sum)
gg$f9010_w <- gg$f9010/gg$sf9010

gg$sf9910 <- ave(gg$f9910,gg$year,FUN=sum)
gg$f9910_w <- gg$f9910/gg$sf9910

gg$lnbwkrs <- ifelse(gg$count>0,log(gg$count),NA)

# gg$leffref <- simplelag(gg$effref,gg$est)

gg$r_lyof <- ifelse(is.na(gg$nblayoffs+gg$effref) | gg$effref==0,0,gg$nblayoffs/gg$effref)
gg$r_lyof <- pmin(gg$r_lyof,1)

gg$r_lyof_cum <- ave(gg$r_lyof,gg$est,FUN=function(x) cumsum(x)) 
gg$lyof_cum <- ave(gg$r_lyof>0,gg$est,FUN=function(x) cumsum(x)) 
# 
# gg$n_lyof_cum <- ave(ifelse(is.na(gg$nblayoffs),0,gg$nblayoffs),gg$est,FUN=function(x) cumsum(x)) 
# summary(gg$n_lyof_cum)

gg$firm <- substr(gg$est,1,9)

gg$sf90102 <- ave(gg$f9010*gg$pds,gg$year,FUN=function(x) sum(x,na.rm=T))
gg$f9010_w2 <- gg$f9010*gg$pds/gg$sf90102


gg$sf99102 <- ave(gg$f9910*gg$pds,gg$year,FUN=function(x) sum(x,na.rm=T))
gg$f9910_w2 <- gg$f9910*gg$pds/gg$sf99102

# gg$ynblayoffs <- gg$year-ave(ifelse((gg$nblayoffs>0)*gg$year>0,(gg$nblayoffs>0)*gg$year,10000),gg$est,FUN=function(x) min(x,na.rm=T))
# gg$miny_nblayoffs <- ave(ifelse((gg$nblayoffs>0)*gg$year>0,(gg$nblayoffs>0)*gg$year,10000),gg$est,FUN=function(x) min(x,na.rm=T))
# gg$sumy_nblayoffs <- ave((gg$nblayoffs>0),gg$est,FUN=function(x) sum(x,na.rm=T))
# gg$nblayoffs2 <- ifelse(gg$miny_nblayoffs==gg$year,0,gg$nblayoffs)
# gg$miny_nblayoffs2 <- ave(ifelse((gg$nblayoffs2>0)*gg$year>0,(gg$nblayoffs2>0)*gg$year,10000),gg$est,FUN=function(x) min(x,na.rm=T))

saveRDS(gg,"gg.rds")


#-----------------------------------------------------------------------------#
#### DATA IMPORT ONCE CREATED  ####
#-----------------------------------------------------------------------------#

rm(gg)
gg <- data.table(readRDS("gg.rds"))

gg <- gg[order(gg$est,gg$year),]
gg[,c("sf9910","f9910_w","f9910xf9910",
"f9099","f9910","sf9010","lndnbwkrs_neg",
"lcount","lnbwkrs"):=NULL]

gg$f9010_min <- ave(gg$f9010,gg$est,FUN=min)
field <- gg$f9010_w2>0 & is.na(gg$f9010_w2)==F & gg$f9010_min>1  & is.na(gg$f9010_min)==F

# table(gg$f9010_min)

#### Descriptives 1 ####
tt0 <- tapply(gg$nblayoffs==0,gg$year,sum)
tt1 <- tapply(gg$nblayoffs>0,gg$year,sum)
tt2 <- tapply(gg$effref,gg$year,sum,na.rm=T)
tt3 <- tapply(gg$nblayoffs,gg$year,sum)
ttp0 <- tapply(gg$pds*(gg$nblayoffs==0),gg$year,sum,na.rm=T)
ttp1 <- tapply(gg$pds*(gg$nblayoffs>0),gg$year,sum,na.rm=T)
ttp2 <- tapply(gg$pds*(gg$effref),gg$year,sum,na.rm=T)
ttp3 <- tapply(gg$pds*(gg$nblayoffs),gg$year,sum,na.rm=T)
ttp4 <- tapply(gg$f9010_w2*gg$f9010xf9010,gg$year,sum,na.rm=T)
ttp5 <- tapply(gg$f9010_w2,gg$year,sum,na.rm=T)

out <- data.frame(tt0,tt1,tt2,tt3,ttp0,ttp1,ttp2,ttp3,ttp4/ttp5)
file.remove("layout_tab.html")
HTML(out,file="layout_tab.html" )

gg <- gg[gg$year<2015,]



gg$nl2002 <- ave((gg$year==2002)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2003 <- ave((gg$year==2003)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2004 <- ave((gg$year==2004)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2005 <- ave((gg$year==2005)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2006 <- ave((gg$year==2006)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2007 <- ave((gg$year==2007)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2008 <- ave((gg$year==2008)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2009 <- ave((gg$year==2009)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2010 <- ave((gg$year==2010)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2011 <- ave((gg$year==2011)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2012 <- ave((gg$year==2012)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2013 <- ave((gg$year==2013)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2014 <- ave((gg$year==2014)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2015 <- ave((gg$year==2015)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2016 <- ave((gg$year==2016)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))
gg$nl2017 <- ave((gg$year==2017)*gg$nblayoffs,gg$est,FUN=function(x) sum(x,na.rm=T))

gg$mynblayoffs_4 <- 1 * (
  ((gg$year-(gg$nl2002>0)*2002)<=-4) + 
  ((gg$year-(gg$nl2003>0)*2003)<=-4) + 
  ((gg$year-(gg$nl2004>0)*2004)<=-4) + 
  ((gg$year-(gg$nl2005>0)*2005)<=-4) + 
  ((gg$year-(gg$nl2006>0)*2006)<=-4) + 
  ((gg$year-(gg$nl2007>0)*2007)<=-4) + 
  ((gg$year-(gg$nl2008>0)*2008)<=-4) + 
  ((gg$year-(gg$nl2009>0)*2009)<=-4) + 
  ((gg$year-(gg$nl2010>0)*2010)<=-4) + 
  ((gg$year-(gg$nl2011>0)*2011)<=-4) + 
  ((gg$year-(gg$nl2012>0)*2012)<=-4) + 
  ((gg$year-(gg$nl2013>0)*2013)<=-4) +
  ((gg$year-(gg$nl2014>0)*2014)<=-4) +
  ((gg$year-(gg$nl2015>0)*2015)<=-4) +
  ((gg$year-(gg$nl2016>0)*2016)<=-4) +
  ((gg$year-(gg$nl2017>0)*2017)<=-4) 
  )

gg$mynblayoffs_3 <- 1 *( 
  ((gg$year-(gg$nl2002>0)*2002)==-3) | 
  ((gg$year-(gg$nl2003>0)*2003)==-3) |     
  ((gg$year-(gg$nl2004>0)*2004)==-3) | 
  ((gg$year-(gg$nl2005>0)*2005)==-3) | 
  ((gg$year-(gg$nl2006>0)*2006)==-3) | 
  ((gg$year-(gg$nl2007>0)*2007)==-3) | 
  ((gg$year-(gg$nl2008>0)*2008)==-3) | 
  ((gg$year-(gg$nl2009>0)*2009)==-3) | 
  ((gg$year-(gg$nl2010>0)*2010)==-3) | 
  ((gg$year-(gg$nl2011>0)*2011)==-3) | 
  ((gg$year-(gg$nl2012>0)*2012)==-3) | 
  ((gg$year-(gg$nl2013>0)*2013)==-3) |
  ((gg$year-(gg$nl2014>0)*2014)==-3) |
  ((gg$year-(gg$nl2015>0)*2015)==-3) |
  ((gg$year-(gg$nl2016>0)*2016)==-3) |
  ((gg$year-(gg$nl2017>0)*2017)==-3))


gg$mynblayoffs_2 <- 1 * (
  ((gg$year-(gg$nl2002>0)*2002)==-2) | 
  ((gg$year-(gg$nl2003>0)*2003)==-2) | 
  ((gg$year-(gg$nl2004>0)*2004)==-2) | 
  ((gg$year-(gg$nl2005>0)*2005)==-2) | 
  ((gg$year-(gg$nl2006>0)*2006)==-2) | 
  ((gg$year-(gg$nl2007>0)*2007)==-2) | 
  ((gg$year-(gg$nl2008>0)*2008)==-2) | 
  ((gg$year-(gg$nl2009>0)*2009)==-2) | 
  ((gg$year-(gg$nl2010>0)*2010)==-2) | 
  ((gg$year-(gg$nl2011>0)*2011)==-2) | 
  ((gg$year-(gg$nl2012>0)*2012)==-2) | 
  ((gg$year-(gg$nl2013>0)*2013)==-2) |
  ((gg$year-(gg$nl2014>0)*2014)==-2) |
  ((gg$year-(gg$nl2015>0)*2015)==-2) |
  ((gg$year-(gg$nl2016>0)*2016)==-2) |
  ((gg$year-(gg$nl2017>0)*2017)==-2)
  )

gg$mynblayoffs_1 <-   1 *(
  ((gg$year-(gg$nl2002>0)*2002)==-1) | 
  ((gg$year-(gg$nl2003>0)*2003)==-1) | 
  ((gg$year-(gg$nl2004>0)*2004)==-1) | 
  ((gg$year-(gg$nl2005>0)*2005)==-1) | 
  ((gg$year-(gg$nl2006>0)*2006)==-1) | 
  ((gg$year-(gg$nl2007>0)*2007)==-1) | 
  ((gg$year-(gg$nl2008>0)*2008)==-1) | 
  ((gg$year-(gg$nl2009>0)*2009)==-1) | 
  ((gg$year-(gg$nl2010>0)*2010)==-1) | 
  ((gg$year-(gg$nl2011>0)*2011)==-1) | 
  ((gg$year-(gg$nl2012>0)*2012)==-1) | 
  ((gg$year-(gg$nl2013>0)*2013)==-1) |
  ((gg$year-(gg$nl2014>0)*2014)==-1) |
  ((gg$year-(gg$nl2015>0)*2015)==-1) |
  ((gg$year-(gg$nl2016>0)*2016)==-1) |
  ((gg$year-(gg$nl2017>0)*2017)==-1)
)

gg$mynblayoffs0 <-   1 * (
  ((gg$year-(gg$nl2002>0)*2002)==0) | 
  ((gg$year-(gg$nl2003>0)*2003)==0) | 
  ((gg$year-(gg$nl2004>0)*2004)==0) | 
  ((gg$year-(gg$nl2005>0)*2005)==0) | 
  ((gg$year-(gg$nl2006>0)*2006)==0) | 
  ((gg$year-(gg$nl2007>0)*2007)==0) | 
  ((gg$year-(gg$nl2008>0)*2008)==0) | 
  ((gg$year-(gg$nl2009>0)*2009)==0) | 
  ((gg$year-(gg$nl2010>0)*2010)==0) | 
  ((gg$year-(gg$nl2011>0)*2011)==0) | 
  ((gg$year-(gg$nl2012>0)*2012)==0) | 
  ((gg$year-(gg$nl2013>0)*2013)==0) |
  ((gg$year-(gg$nl2014>0)*2014)==0) |
  ((gg$year-(gg$nl2015>0)*2015)==0) |
  ((gg$year-(gg$nl2016>0)*2016)==0) |
  ((gg$year-(gg$nl2017>0)*2017)==0)
)

gg$mynblayoffs1 <-   1 * (
  ((gg$year-(gg$nl2002>0)*2002)==1) | 
  ((gg$year-(gg$nl2003>0)*2003)==1) | 
  ((gg$year-(gg$nl2004>0)*2004)==1) | 
  ((gg$year-(gg$nl2005>0)*2005)==1) | 
  ((gg$year-(gg$nl2006>0)*2006)==1) | 
  ((gg$year-(gg$nl2007>0)*2007)==1) | 
  ((gg$year-(gg$nl2008>0)*2008)==1) | 
  ((gg$year-(gg$nl2009>0)*2009)==1) | 
  ((gg$year-(gg$nl2010>0)*2010)==1) | 
  ((gg$year-(gg$nl2011>0)*2011)==1) | 
  ((gg$year-(gg$nl2012>0)*2012)==1) | 
  ((gg$year-(gg$nl2013>0)*2013)==1) |
  ((gg$year-(gg$nl2014>0)*2014)==1) |
  ((gg$year-(gg$nl2015>0)*2015)==1) |
  ((gg$year-(gg$nl2016>0)*2016)==1) |
  ((gg$year-(gg$nl2017>0)*2017)==1)
)

gg$mynblayoffs2 <-   1 * (
  ((gg$year-(gg$nl2002>0)*2002)==2) | 
  ((gg$year-(gg$nl2003>0)*2003)==2) | 
  ((gg$year-(gg$nl2004>0)*2004)==2) | 
  ((gg$year-(gg$nl2005>0)*2005)==2) | 
  ((gg$year-(gg$nl2006>0)*2006)==2) | 
  ((gg$year-(gg$nl2007>0)*2007)==2) | 
  ((gg$year-(gg$nl2008>0)*2008)==2) | 
  ((gg$year-(gg$nl2009>0)*2009)==2) | 
  ((gg$year-(gg$nl2010>0)*2010)==2) | 
  ((gg$year-(gg$nl2011>0)*2011)==2) | 
  ((gg$year-(gg$nl2012>0)*2012)==2) | 
  ((gg$year-(gg$nl2013>0)*2013)==2) |
  ((gg$year-(gg$nl2014>0)*2014)==2) |
  ((gg$year-(gg$nl2015>0)*2015)==2) |
  ((gg$year-(gg$nl2016>0)*2016)==2) |
  ((gg$year-(gg$nl2017>0)*2017)==2))

gg$mynblayoffs3 <- 1 * (
  ((gg$year-(gg$nl2002>0)*2002)==3) | 
  ((gg$year-(gg$nl2003>0)*2003)==3) | 
  ((gg$year-(gg$nl2004>0)*2004)==3) | 
  ((gg$year-(gg$nl2005>0)*2005)==3) | 
  ((gg$year-(gg$nl2006>0)*2006)==3) | 
  ((gg$year-(gg$nl2007>0)*2007)==3) | 
  ((gg$year-(gg$nl2008>0)*2008)==3) | 
  ((gg$year-(gg$nl2009>0)*2009)==3) | 
  ((gg$year-(gg$nl2010>0)*2010)==3) | 
  ((gg$year-(gg$nl2011>0)*2011)==3) | 
  ((gg$year-(gg$nl2012>0)*2012)==3) | 
  ((gg$year-(gg$nl2013>0)*2013)==3) |
  ((gg$year-(gg$nl2014>0)*2014)==3) |
  ((gg$year-(gg$nl2015>0)*2015)==3) |
  ((gg$year-(gg$nl2016>0)*2016)==3) |
  ((gg$year-(gg$nl2017>0)*2017)==3)
  )

gg$mynblayoffs4 <- 1* (
  ((gg$year-(gg$nl2002>0)*2002) %in% c(4:20)) + 
  ((gg$year-(gg$nl2003>0)*2003) %in% c(4:20)) +
  ((gg$year-(gg$nl2004>0)*2004) %in% c(4:20)) +
  ((gg$year-(gg$nl2005>0)*2005) %in% c(4:20)) +
  ((gg$year-(gg$nl2006>0)*2006) %in% c(4:20)) +
  ((gg$year-(gg$nl2007>0)*2007) %in% c(4:20)) +
  ((gg$year-(gg$nl2008>0)*2008) %in% c(4:20)) +
  ((gg$year-(gg$nl2009>0)*2009) %in% c(4:20)) +
  ((gg$year-(gg$nl2010>0)*2010) %in% c(4:20)) +
  ((gg$year-(gg$nl2011>0)*2011) %in% c(4:20)) +
  ((gg$year-(gg$nl2012>0)*2012) %in% c(4:20)) +
  ((gg$year-(gg$nl2013>0)*2013) %in% c(4:20)) +
  ((gg$year-(gg$nl2014>0)*2014) %in% c(4:20)) +
  ((gg$year-(gg$nl2015>0)*2015) %in% c(4:20)) +
  ((gg$year-(gg$nl2016>0)*2016) %in% c(4:20)) +
  ((gg$year-(gg$nl2017>0)*2017) %in% c(4:20))
)


gg[,c("nl2002","nl2003","nl2004","nl2005","nl2006","nl2007","nl2008","nl2009","nl2010","nl2011",
      "nl2012","nl2013","nl2014","nl2015","nl2016","nl2017"):=NULL]

str(gg)
gg$lyof <- (gg$r_lyof>0)*1

#### descriptives 2 ####
field <- gg$f9010_w2>0 & is.na(gg$f9010_w2)==F & gg$f9010_min>1  & is.na(gg$f9010_min)==F
gg$lnbwkrs<- log(gg$count)
d1 <- wtd.summary(gg$f9010xf9010[field],w=gg$f9010_w2[field])
d2 <- wtd.summary(gg$lnbwkrs[field],w=gg$f9010_w2[field])
d3 <- wtd.summary(gg$lnnbwkrs_cumneg[field],w=gg$f9010_w2[field])
d4 <- wtd.summary(gg$lyof[field],w=gg$f9010_w2[field])
# d4b <- wtd.summary(gg$lyof,w=(gg$pds*gg$count))
# d4c <- wtd.summary(gg$lyof[gg$f9010_w2>0],w=(gg$pds[gg$f9010_w2>0]))
# d4d <- wtd.summary(gg$lyof[gg$f9010_w2>0],w=(gg$pds[gg$f9010_w2>0]*gg$count[gg$f9010_w2>0]))
d5 <- wtd.summary(gg$r_lyof[field],w=gg$f9010_w2[field])
d6 <- wtd.summary(gg$r_lyof_cum[field],w=gg$f9010_w2[field])
d7 <- wtd.summary(gg$lyof_cum[field],w=gg$f9010_w2[field])

des_lyof <- rbind(d1,d2,d3,d4,d5,d6,d7) 
rownames(des_lyof) <- NULL

des_lyof
write.csv(des_lyof,"des_lyof_2.csv")


#-----------------------------------------------------------------------------#
#### simple cum layoffs events models  ####
#-----------------------------------------------------------------------------#

# Size restriction
field <- gg$f9010_w2>0 & is.na(gg$f9010_w2)==F & gg$f9010_min>1 & is.na(gg$f9010_min)==F
gg <- gg[field==T,]
rownames(gg) <- NULL
gc()

mod90_1 <- felm(I(100*f9010xf9010) ~ 
                  factor(year)
                + lyof_cum
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg)
summary(mod90_1)
screenreg(mod90_1)

mod90_2 <- felm(I(100*f9010xf9010) ~ 
                  factor(year)
                + lyof_cum
                + log(count)
                + lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg )
summary(mod90_2)

screenreg(list(mod90_1,mod90_2),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_1,mod90_2),stars=c(0.1,0.05,0.01),file="layoff_cum.html",digits=3)

gc()
rm(mod90_1,mod90_2)

#-----------------------------------------------------------------------------#
#### Cumulative layoffs models  ####
#-----------------------------------------------------------------------------#

mod90_3 <- felm(I(100*f9010xf9010) ~ 
                  factor(year)
                + r_lyof_cum
                |est|0|est, 
                weights=gg$f9010_w2, 
                data=gg)
summary(mod90_3)
screenreg(mod90_3)

mod90_4 <- felm(I(100*f9010xf9010) ~ 
                  factor(year)
                + r_lyof_cum
                + log(count)
                + lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg )
summary(mod90_4)

screenreg(list(mod90_3,mod90_4),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_3,mod90_4),stars=c(0.1,0.05,0.01),file="r_layoff_cum.html",digits=3)

rm(mod90_3,mod90_4)


#-----------------------------------------------------------------------------#
#### MULTI EVENTS MODELS  ####
#-----------------------------------------------------------------------------#

mod90_5 <- felm(I(100*f9010xf9010) ~ 
                  factor(year)+
                #   + mynblayoffs_6  
                # +mynblayoffs_5
                +mynblayoffs_4
                +mynblayoffs_3
                +mynblayoffs_2
                # +mynblayoffs_1
                +mynblayoffs0
                +mynblayoffs1
                +mynblayoffs2
                +mynblayoffs3
                +mynblayoffs4
                # +mynblayoffs5
                # +mynblayoffs6
                # +lnbwkrs
                # +lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg)
summary(mod90_5)
screenreg(mod90_5)

mod90_6 <- felm(I(100*f9010xf9010) ~ 
                  factor(year)+
                #   + mynblayoffs_6  
                # +mynblayoffs_5
                +mynblayoffs_4
                +mynblayoffs_3
                +mynblayoffs_2
                # +mynblayoffs_1
                +mynblayoffs0
                +mynblayoffs1
                +mynblayoffs2
                +mynblayoffs3
                +mynblayoffs4
                # +mynblayoffs5
                # +mynblayoffs6
                +log(count)
                +lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg )
summary(mod90_6)

screenreg(list(mod90_5,mod90_6),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_5,mod90_6),stars=c(0.1,0.05,0.01),file="did_mult_layoff_2.html",digits=3)

rm(mod90_5,mod90_6)
#-----------------------------------------------------------------------------#
#### Trends  ####
#-----------------------------------------------------------------------------#

mod90_7 <- felm(I(100*f9010xf9010) ~ 
                  year
                #   #   + mynblayoffs_6  
                #   # +mynblayoffs_5
                #   +mynblayoffs_4
                # +mynblayoffs_3
                # +mynblayoffs_2
                # # +mynblayoffs_1
                # +mynblayoffs0
                # +mynblayoffs1
                # +mynblayoffs2
                # +mynblayoffs3
                # +mynblayoffs4
                # # +mynblayoffs5
                # # +mynblayoffs6
                # # +lnbwkrs
                # # +lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg)
summary(mod90_7)


mod90_8 <- felm(I(100*f9010xf9010) ~ 
                  year+
                  #   + mynblayoffs_6  
                  # +mynblayoffs_5
                  +mynblayoffs_4
                +mynblayoffs_3
                +mynblayoffs_2
                # +mynblayoffs_1
                +mynblayoffs0
                +mynblayoffs1
                +mynblayoffs2
                +mynblayoffs3
                +mynblayoffs4
                # +mynblayoffs5
                # +mynblayoffs6
                # +lnbwkrs
                # +lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg)
summary(mod90_8)

mod90_9 <- felm(I(100*f9010xf9010) ~ 
                  year+
                  r_lyof_cum
                  #   + mynblayoffs_6  
                  # +mynblayoffs_5
                #   +mynblayoffs_4
                # +mynblayoffs_3
                # +mynblayoffs_2
                # # +mynblayoffs_1
                # +mynblayoffs0
                # +mynblayoffs1
                # +mynblayoffs2
                # +mynblayoffs3
                # +mynblayoffs4
                # # +mynblayoffs5
                # # +mynblayoffs6
                # # +lnbwkrs
                # # +lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg)
summary(mod90_9)
mod90_10 <- felm(I(100*f9010xf9010) ~ 
                  year+
                  lyof_cum
                #   + mynblayoffs_6  
                # +mynblayoffs_5
                #   +mynblayoffs_4
                # +mynblayoffs_3
                # +mynblayoffs_2
                # # +mynblayoffs_1
                # +mynblayoffs0
                # +mynblayoffs1
                # +mynblayoffs2
                # +mynblayoffs3
                # +mynblayoffs4
                # # +mynblayoffs5
                # # +mynblayoffs6
                # # +lnbwkrs
                # # +lnnbwkrs_cumneg
                |est|0|firm, 
                weights=gg$f9010_w2, 
                data=gg)
summary(mod90_10)

screenreg(list(mod90_7,mod90_8,mod90_9,mod90_10),stars=c(0.1,0.05,0.01),digits=3)
htmlreg(list(mod90_7,mod90_8,mod90_9,mod90_10),stars=c(0.1,0.05,0.01),digits=3,
        file="did_mult_layoff_2_trend.html")

rm(mod90_7,mod90_8,mod90_9,mod90_10)